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'S ■ ABSTRACT 

I We study the reaction of a globular star cluster to a time-varying tidal perturbation 

^ ■ (gravitational shock) using self-consistent N-body simulations and address two questions. 



> 






First, to what extent is the cluster interior protected by adiabatic invariants. Second, 
how much further energy change does the postshock evolution of the cluster potential 



SO ! produce and how much does it affect the dispersion of stellar energies. We introduce 

^ ■ the adiabatic correction as ratio of the energy change, {/^E) , to its value in the impulse 

^ . approximation. When the potential is kept fixed, the numerical results for the adiabatic 

O ! correction for stars with orbital frequency u) can be approximated as (1 -(-a;^r^)~'^. For 



shocks with the characteristic duration of the order the half-mass dynamical time of the 

cluster, r ^ tdyn,h, the exponent 7 = 5/2. For more prolonged shocks, r ^ 4:tdyn,h, the 

p . adiabatic correction is shallower, 7 = 3/2. When we allow for self-gravity and potential 

O 



oscillations which follow the shock, the energy of stars in the core changes significantly, 

Jh I while the total energy of the system is conserved. Paradoxically, the postshock potential 

^ ■ fluctuations reduce the total amount of energy dispersion, (AE"^). The effect is small but 

real and is due to the postshock energy change being statistically anti-correlated with 

the shock induced heating. These results are to be applied to Fokker-Planck models of 

/\ ' the evolution of globular clusters. 

Subject headings: stellar dynamics — globular clusters: general — methods: numerical 



1. Introduction 

The motion of star clusters in the Galaxy is pri- 
marily determined by the monopole component of the 
Galactic gravitational force; their internal dynamics 
is often influenced by the tidal force. A steady tidal 
force imposes a cutoff of the stellar distribution in the 
cluster at the tidal radius, whereas fast encounters 
with th e Galactic disk or bulge cause gra vitational 
shocks ( Ostriker, Spitzer, fc Chevalier 1972 , hereafter 
OSC; ISpitzer 1987| ). 

The overall effect of tidal shocks is the enhanced 
evaporation of clusters via the following processes. 
First, stars on average gain energy {{AE) > 0) reduc- 
ing the cluster binding energy. After a short period of 
contraction immediately following the shock, the clus- 
ter as a whole expands and some stars find themselves 
outside the tidal boundary. They are now captured by 
the galactic potential and are effectively lost from the 
cluster. Second, tidal shocks induce a dispersion of 
stellar energies {{AE'^) > 0), named the "tidal shock 
relaxation" , which is very similar to the diffusion in 



energy space due to two-body relaxation (Kundic & 
Ostriker 1995; hereafter KO). Both of these processes 
supply stars to the high-energy tail of the velocity 
distribution, which is not bound by the cluster po- 
tential. Those stars eventually escape, thus enhancing 
the evaporation of the cluster. Finally, shock-induced 
relaxation speeds up core collapse, which in turn leads 
to an ever faster evolution, faster evaporation, and 
subsequent disintegration of the cluster. 

OSC considered the problem in the impulse ap- 
proximation, wherein the stars are assumed to move 
little during the shock. This is appropriate in the 
outer regions of the cluster, near the tidal radius. In 
the cluster core, the impact of shocks is very differ- 
ent. The rapidly orbiting stars smear the perturba- 
tion over many rotations around the center of the clus- 
ter. Because the conservation of adiabatic invariants 
prevents stars from gaining energy, the shocking ef- 
fect is significantly reduced. Spitzer (1987) studied 
this problem in the harmonic potential approxima- 
tion and found that (AE) is exponentially suppressed. 
Recently Weinberg (1994a) suggested that the actual 
impact of the shock is stronger because of the numer- 
ous resonances in a system with more than one degree 
of freedom. Since the original formula for (AE) was 
derived by OSC in the impulse approximation, we de- 
fine the adiabatic correction as the ratio of the actual 
energy change to its "impulsive" value. 



In addition to the remaining question concern- 
ing the proper treatment of the adiabatic correction, 
there is a much larger open question. After the per- 
turbing force has ceased, the cluster finds itself out 
of equilibrium. Then there will be an interval of sev- 
eral dynamical times while the cluster oscillates until 
phase mixing and Landau damping bring it into a 
new equilibrium state. During this phase, the total 
energy of the cluster will be conserved but individ- 
ual stellar orbits may gain or lose energy in response 
to the fiuctuating gravitational field. This may pro- 
duce additional both the first and second order energy 
changes, {AE) and {AE-^). These secondary effects 
have been studied very little and may contribute to 
the cluster evolution, along with the better known 
primary effects. 

As an instructive exercise, consider how the stars 
gain energy at the expense of the binding energy of the 
cluster in an impulsive shock. Let Ei be the energy 
of the ith star, 
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where mj, Fj, and v^ are the mass, position, and ve- 
locity of the star, respectively. The total energy of 
the system is a sum of the kinetic (T) and potential 
(W) terms: 
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This can be rewritten as 
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Immediately after the shock, the stars have gained 
some kinetic energy, AEi^sh, but have not yet had 
time to change the potential. As the result, the energy 
of the cluster changes by 



AEt. 



E^i^. 
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After the period of virialization following the shock, 
both the kinetic and the potential energies of the clus- 
ter change, AWvir — —2 ATVir, but the total energy is 
conserved. This implies AWvir = 2AEtot, and com- 
bined with equations (||) and (Q) gives the final energy 
change of the stars: 

Y, A^^,vir = AEtot + AWvir = ^ 3 A^,,,h. (5) 



The stars gain on average three times the initial kick 
due to the shock. The extra energy comes from the 
reduction of the binding energy of the cluster. A 
schematic picture of the potential change is presented 
in Figure U 

We use self-consistent N-body simulations to in- 
vestigate the adiabatic correction as a function of the 
star energy and position in the cluster. This tech- 
nique allows us to study directly the self-gravitating 
response of the cluster to the tidal perturbation, ne- 
glected in prior models. We find that collective effects 
play an important role in transporting the pertur- 
bation from the outer layers deep inside the cluster. 
While the postshock evolution conserves the total en- 
ergy, the energy of individual stars changes because 
of the readjustment of the potential. 

A previous approach to include tidal shocks in the 
evolutionary Fokker-Planck code ( [Gnedin fc Ostriker 
1997) has been to first input the energy changes asso- 
ciated with the shock, and then update the potential 
using the adiabatic invariants. These adiabatic invari- 
ants change as a result of shocking and so the above 
procedure is not exactly correct, but no better alter- 
native was available prior to the current work. Also, 
this method does not include the energy transfer be- 
tween stars subsequent to the shock which gives rise 
to additional changes in {AE) and {AE"^). Our cur- 
rent approach is fully self-consistent and treats all of 
the above processes exactly. 

A direct application of the linear perturbation the- 
ory has been done by Weinberg (1994b, c) and Murah 
& ' Veinberg (1997a-c). Both methods, our fully nu- 
merical and the semi-analytical by the authors, lead 
to the correct results but differ in implementation. 

We review the analytic calculations of the adiabatic 
correction of Spitzer and Weinberg in §0. In §0, we 
describe the N-body code and the parameters of the 
star cluster. In §0, we present the results of our sim- 
ulations of disk shocking and compare them with the 
analytical estimates. We consider the problem twice, 
keeping the cluster potential fixed, and allowing for 
self-gravity. Also in §0, we investigate radial shocking 
to test the effects of the self-gravitating response of 
the cluster noted above. We summarize and discuss 
the implications of our results for the Fokker-Planck 
calculations of the dynamical evolution of stellar sys- 
tems in §|6[ 



2. Adiabatic Corrections 

We need a parameter to quantify the transition 
from the impulse to the adiabatic regime of tidal 
shocking. Let r be the characteristic duration of the 
shock. We define the adiabatic parameter 



(6) 



where lo is the stellar orbital frequency, or, defined 
in terms of the root-mean-square velocity Vrms of 
stars at a distance r from the center of the cluster, 
uj{r) = Vrms{r)/r. This frequency is directly related 
to the (inverse) dynamical time of stars at a given 
distance r. The effective time to cross the disk of the 
Galaxy of half-thickness H for the cluster with the 
perpendicular velocity V is 
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For the shock duration r representative for globular 
clusters in the Galaxy, x ^ 1 close to the tidal radius, 
and a; ^ 1 in the core. To the second order in per- 
turbation, the initial energy change inside the cluster 
relative to the outer parts is independent of the shock 
amplitude. Therefore we can define a function A{x) 
of a single variable x as the ratio of the actual energy 
change to its impulse value. The adiabatic correction 
A{x) should satisfy the following criteria: 

A{x) = 1, for x < 1; A{x) < 1, for x > 1. (8) 

In the harmonic approximation, the first and sec- 
ond order terms for the average energy change of stars 
at the distance z from the equatorial plane of the clus- 
ter are given by Spitzer (1987) and KO, respectively: 



(Ai?) = i((At;,)2) = 
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where g^ is the maximum vertical gravitational ac- 
celeration produced by the disk, and (z^) = r^/S for 
a phase-mixed spherical mass distribution. Here 



Ai{x) = exp(-2x^) , 
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A2{x) = - exp(-2x2) 
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and we denote by the subscripts 1 and 2 the adia- 
batic correction for the first and second order terms. 



respectively. If we put Ai = A2 = I, we recover the 
impulse approximation result. Note that the value of 
A2 jumps by the factor 9/5 as a; ^ 0; i.e., the har- 
monic approximation enhances the shock impact at 
large radii. However, we should not expect it to be 
unity, since the assumption of the parabolic potential 
is not valid on the cluster periphery. 

The derivation of the above result assumes that 
the oscillation frequencies of the system are not com- 
mensurable with the perturbation frequency. In gen- 
eral, it has been proved that under this condition the 
adiabatic factor vanishes faster than any power of x 



and finally obtain 



(Kruskal 1962), consistent with the exponential vari- 
ation (Spitzer 1987, p. 



48). 



Recently Weinberg (1994a) showed that resonances 
do typically occur in a system with more than one de- 
gree of freedom. If the system is represented by a com- 
bination of multi-dimensional nonlinear oscillators, it 
becomes likely that some of the perturbation frequen- 
cies will be commensurable with the low oscillation 
frequencies of the stars. Then those stars receive a 
significant kick from the perturbation and no longer 
conserve their actions. Averaging over an ensemble 
of stars can give an appreciable energy change. As a 
result of summing over the resonant terms, the adi- 
abatic factor A{x) is not exponentially small for the 
large values of x, but rather follows approximately a 
power-law form. 

Formulae derived by Weinberg (1994b; for exam- 
ple, his eq. [24]) contain multiple infinite sums and 
integrals over all frequency range allowed to the sys- 
tem. Instead of a direct application of the linear the- 
ory formalism, we derive an asymptotic result for the 
long duration of the shock, r, and fit the intermedi- 
ate reginie^_We follow closely the suggestion by ^ 



Trc maine (1996)[ Consider a shock with the dura- 
tion comparable to or larger than the half-mass ra- 
dius dynamical time of the cluster. The acceleration 
of stars due to the perturbation is roughly g^ z/H, 
and their energies change significantly only over the 
period when the stars are in resonance with the per- 
turbation, Ai ^ uj~^. Thus those stars would gain 
about 
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The energy change per unit mass is A£' — AEres'x 
(fraction of resonant stars). Equations (19) and (24) 
of Weinberg (1994b) indicate that the number of stars 
at the peak amplitude scales as 1/r. Therefore, we 
approximate the fraction of resonant stars as (w t)~^ 
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where we have substituted equation (|7|) for r. It fol- 
lows from the equation above that the asymptotic 
form of the adiabatic correction is 



A{x) oc X 



for a: ::^ 1. 



(13) 



Now we can construct the correction factor for all val- 
ues of X from its asymptotic behavior. The simplest 
fitting formula is 



A{x)^ {l + x') 



-3/2 



(14) 



In the following, we call equation (|lj) the Wein- 
berg correction, as it was derived from considerations 
brou ght f orward by M. Weinberg.^] We refer to equa- 
tion ( |lOa| ) as the Spitzer correction. We plot the two 
adiabatic corrections on Figure 0. The exponential 
drops much faster than the power-law in the central 
region of the cluster. 

3. Self-Consistent Field Code 

We perform N-body simulations of a single tidal 
shock in a spherical system in order to investigate 
the validity of the adiabatic corrections considered in 
the previous section. We use first the self-consistent 



field (SCF) method described in Hernquist & Ostriker 
(1992). The code is designed to reduce the numerical 
relaxation, which arises from close two-body interac- 
tions. Instead of direct calculation of forces for all 
particles, the code computes orbits of the particles in 
a smooth potential of the cluster, expanded in a se- 
ries of predetermined basis functions. The coefficients 
of the expansion are updated each time step based 
on the particle positions. This procedure provides a 
self-consistent solution to the Poisson equation. The 
details of th e calculation of the ex pansion coefhcients 
are given in Hernquist fc Ostrikei (1992), whose set of 
basis functions we use here. We have done simulations 
with two sets of the expansion coefficients (rimax = 6, 
Imax = 4) and {rimax = 10, Imax = 6). In both cases 
the results are indistinguishable and therefore are in- 
dependent of the numerical method. 



^ Note that equation (hj) is an asymptotic fitting formula for 
the case uit "^ 1 . T he linear tlieory of Weinberg can be used 
( Murali & Weinberg 1997a-c) in all regimes, albeit the resulting 
expressions are rather complex. 



As our initial model for the cluster we take the 
King model ( pK^ing 1966 ) with the structural param- 
eter Wq — 4 (corresponding to the concentration 
c = 0.84). This is a relatively loosely bound clus- 
ter, for which the effects of tidal shocks should be 
important. Working units of the code are such that 
G = M = Re ~ I, where M is the total mass of the 
cluster and Re is the core radius. The tidal radius is 
determined by the concentration; Rt « 6.92 i?c- The 
characteristic time for the whole cluster is the dy- 
na mical time a t the half-mass radius Rh ( Binney fc 
Trc jmainc 19"87| , eq. [2-30]): 
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in the code units. For our cluster, R/^ w 1.58i?c- We 
have chosen a time step which is 1% of the half-mass 
dynamical time, At = 0.01 tdyn,h- This time step is 
sufficient to calculate accurately stellar orbits even at 
the center, since the shortest orbital period is Pc = 

^tdyn,c = '^{ph/pcY''^tdyn,h ~ 1.26tdyn.h, where Pc is 
the central density, and ph is the density at half-mass 
radius; ph ~ O.lpc for this cluster. 

We use 10^ equal-mass particles to model the clus- 
ter. The number of stars in observed globular clusters 
is of the same order, so that effects of "noise" (such 
as the Poisson noise) are equally present in both our 
simulations and real clusters. The numerical accuracy 
of the calculations is secured by the conservation of 
the total energy of the cluster in isolation at the level 
AE/Er^ 10"^ 

4. Numerical Results 

We investigate first the impulsive shocking that oc- 
curs during one time step. It allows us to study the 
self-gravitating response of the cluster separately from 
the adiabatic corrections. In run A the potential is 
fixed, and in run B it is fully self-consistent. 

Then we allow for the time-varying perturbation 
with the duration of the order of the dynamical time 
of the cluster. For the real clusters in the Galaxy, a 
typical shock lasts r « 1 — Stdynji- Run C explores 
the case with the fixed potential, run D is the final 
self-consistent simulation. The parameters of all runs 
are summarized in Table nl 



4.1. Impulsive perturbation 

Runs A and B study disk shocking in the impulsive 
regime. We apply an impulse of the form 



dt 
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during one time step at the time t — tdyn,h and then 
follow the evolution of the cluster until t = 20 tdynji- 
We take hmp = 1 h^ the code units, which corre- 
sponds to limp ~ 0-2^i~"^^rfJ„ h- Other components 
of the stellar velocity remain unchanged during the 
impulse. The overall effect of the shock is character- 
ized by the relative reduction in the binding energy, 
AEtot/Eund = 0.002. 

Figure S illustrates the time evolution of the sys- 
tem. We plot here the change of the total energy 
of the cluster, as well as the partial contributions of 
the kinetic and potential energies. The total energy 
is very well conserved after the impulse, but both T 
and W undergo several oscillations. The lower panel 
of Figure |3| shows deviations from the virial equilib- 
rium as indicated by the ratios —2 T/W and —2 T/V, 
where V is the Virial of the system 



^-E 



Xfc -Xfc, 



(17) 



and V — ~2T for any system of particles moving in a 
finite region of space with finite velocities, wh en aver- 
aged over time (e.g.. Landau fc Lifshitz 1988 ). Thus, 
for an equilibrium system the virial ratios should be 
unity. 

From these plots we see how the cluster responses 
to the perturbation. At first, it is compressed since all 
stars undergo instant acceleration toward the center. 
Their kinetic energy is on average increased and soon 
they move outward, and an expansion follows the rel- 
atively short period of contraction. Now the stars lose 
kinetic energy, and after about itdyn,h the virial re- 
lations are restored. However, the cluster continues 
to expand following inertia and undergoes a few more 
oscillations. The virial equilibrium is approached in 
about 10— 15 dynamical times. As expected, we find 
ATvir = -AEtot and AW^ir = 2 AEtot to be within 
5% accuracy at the end of the simulation. 

Now we look at the distribution of the energy 
changes inside the cluster. First, we fix the cluster 
potential (run A) and see how well our results agree 
with the analytical predictions. Figure shows the 
relative changes {AE)/E and {AE^)/E^ at the end 



of the simulation as a function of the initial energy, 
Ei. The stars are grouped into 100 bins, each having 
10^ particles. This assures proper averaging over the 
phase space of stellar orbits. We expect that statisti- 
cal errors do not exceed the 1% level. The solid line 
shows the expected energy changes for stars of the 
same initial energy: 
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(18a) 



(Ai?^> = ^""'^f'^^ ^^L.(l + X..),(18b) 

where Ai is the time step, r and Vrms arc the mean 
position and rms velocity of stars of energy E, and 
Xr V is the position- velocity correlation function (see 
§M). The data points lie right along the predicted 
curve, confirming that the binning procedure is cor- 
rect and numerical errors are not significant. 

Next, we study the self-gravitating response of the 
cluster to the perturbation. Figure ra shows the en- 
ergy changes for run B. Right after the shock, both 
the first and second order energy changes follow their 
estimates for the fixed potential. Twenty dynamical 
times later, AE is significantly enhanced everywhere 
in the cluster. The energy dispersion is also much 
higher in the cluster core. What causes such a differ- 
ence? 

First we check whether it is simply a result of 



nu merical re laxation. Recent studies ( Hernquist fc 
Ba: nes 1990 ; Hernquist & Ostriker 1992; Weinberg 
19!! 3|) showed that the self-consistent field method 
does not totally avoid relaxation because of the fi- 
nite number of particles used in calculation of the 
potential. Although reduced, compared to the di- 
rect summation methods, numerical relaxation drives 
dispersion of stellar energies over the time scale pro- 
portional to the number of particles in the system. 
We performed test simulations varying the number of 
particles from A^ = 10^ to A^ = 2 x 10^. The energy 
dispersion, {AE^)n, does decrease with increasing A^, 
indicating that there exists a spurious numerical re- 
laxation. In the core, the results can be approximated 

by 



{AE')n = 
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(19) 



The first order energy change is not affected by 
the numerical effects as the dispersion they induce 
has zero mean. The observed heating is the result 
of the cluster potential attaining a new equilibrium. 



causing energy exchange and phase mixing of stellar 
orbits. 

As a result of the heating, the cluster expands over- 
all. Figure o shows the density distribution at the end 
of run B. Some stars move beyond the tidal radius but 
are not necessarily lost: they are still gravitationally 
bound to the cluster. 

In the new equilibrium, the depth of the potential 
well is decreased relative to the initial configuration 
(cf. Figure |l|). Since the amplitude of the pertur- 
bation is small, the potential changes self-similarly, 
i.e. A$ oc $i, where $i is the initial potential. We 
find that the simulation results can be fitted by an 
almost independent sum of the impulsive heating (as 
calculated by OSC) and the potential gain: 

(A^) - (A^),h + AEpou (20) 

where (A£')sh is given by equation (p^, and 

AEpot = c Il,j, (-$,;), (21) 

where c is a normalization constant such that the sum 
of total energy change over all particles is twice the 
total energy change of the cluster (required by the 
Virial theorem; cf. eq. j|]). Note that AEpot is posi- 
tive for all stars and is unaffected by the phase space 
averaging. 

The evolution of the energy dispersion is less ob- 
vious. Figure o shows that {AE'^) decreases relative 
to the impulsive value in the middle of the cluster. 
We found this behavior generic for all simulations we 
performed. The distribution of energies is illustrated 
in Figure 0, which shows the detailed structure of one 
bin at the half-mass radius. The energy change AE 
is slightly asymmetric with respect to the center of 
the bin, (AE), immediately after the shock but it is 
very well phase mixed at the end of the run. The 
distribution is adequately represented by the Gaus- 
sian form although the dispersion is indeed reduced 
by about 20% relative to the initial impulse. The 
normal shape of the distribution is intrinsic as it is 
unaffected by varying the size of the bins. The skew- 
ness and kurtosis of the distribution are small and 
statistically insignificant. 

The reduction in the dispersion of stellar energies 
after potential fluctuations, or the "cooling" of the 
cluster, can be explained by the anti-correlation of the 
initial energy change and the subsequent change due 
to self-consistent oscillations. In a compressive shock, 
stars that were moving inwards before the shock gain 



energy and those moving outwards lose energy. But, 
after the shock there is a compressive wave moving 
inwards (which subsequently reflects and moves out- 
wards), so those particles which had gained energy 
are in phase with the wave and those which had lost 
energy are out of phase. In an expanding shock, all 
directions are reversed but the relation between the 
energy gain and the similarity of phase is maintained. 
This process induces an anti-correlation between the 
initial energy change and the subsequent change when 
the wave reverses. 

Figure H shows the correlation of the additional 
energy change following the perturbation {AEpat = 
AEf — AEsh) with the radial velocity of stars in the 
half-mass bin. Figure demonstrates that this pro- 
duces an anti-correlation of AEpot with the initial 
kick, AEkin- The phase mixing of stars during po- 
tential fluctuations leads to a very smooth flnal dis- 
tribution of energies in the bin (Figure |7|). Since it 
fits the expected Gaussian function, the distribution 
can be described by only the first two moments, (AE) 
and (AE^). 

Another point of view has been suggested by the 
referee, M. Weinberg (also, Johnston et al. 



Stellar actions, the radial action Ir and the angu- 
lar momentum J, change after the impulsive shock. 
During a relatively slow virialization process follow- 
ing the shock, the actions are conserved and so is the 
phase space available to stars of the same initial en- 
ergy. However, the energy space has decreased after 
the virialization, causing the corresponding decrease 
of the energy dispersion, {AE"^). 

We have checked the postshock conservation of 
the phase space in run B. The angular momentum 
is virtually conserved, although (AJ^) increases in 
the cluster core due to numerical relaxation effects. 
The radial action grows during the virialization, pre- 
sumably because of the change of the potential. The 
dispersion of the radial action, (A/^), increases in the 
core due to numerical effects but decreases slightly in 
the middle parts of the cluster. We find that the value 
of (A/^) drops by 9% at the half-mass radius and by 
13% at larger radn. At the tidal radius, (AJ^) again 
increases relative to the impulsive value. Although 
of lower amplitude than the corresponding reduction 
of the energy dispersion (39% and 34%, respectively), 
the decrease of the dispersion of the radial action dur- 
ing the postshock phase is not a numerical effect. This 
must be a contribution of the collective oscillatory ef- 
fects, in addition to the change of the potential. 



4.2. Time-varying perturbation 

Having explored the general features of the self- 
gravitating response of the cluster to the impulsive 
perturbation, we now turn to modeling the more re- 
alistic time-varying shocking. For the disk shocking 
in the Galaxy, the vertical gravitational acceleration 
towards the Galactic plane can be approximated by 
a Gaussian function of the height Z, with the charac- 
teristic thickness H. If we assume a constant vertical 
component of the cluster velocity during the passage 
through the disk, the tidal force will vary as e~* . 
Therefore, we apply a shock of the form 



dv 
~dt 
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where lexp = hmp At/(y^ t) to give the same total 
energy change for r — )■ as in the impulsive shocking. 
Here r is the effective duration of the shock (eq. 0), 
which we have chosen to be equal to the half-mass dy- 
namical time, T = tdyn,h, for runs C and D. We start 
the simulation long before the maximum amplitude 
of the shock occurs at to = 4 tdyn.h and then follow 
the evolution for iQtdyn.h- 

Figure Iffl shows the total energy change for run 
D. This plot is similar to Figure for the impulsive 
shocking, except that here AEtot is reduced because 
of the conservation of the adiabatic invariants of cen- 



tral stars (see also §4.3) 



The distribution of the energy changes for the fixed 
potential (Figure |l^) is intermediate of those pre- 
dicted by Spitzer's (eq. |lO|) and Weinberg's (eq. 
iQ) adiabatic corrections. We find that the results 
in a new equilibrium can be fitted with a reasonable 
accuracy by equations (|l^) with the adiabatic correc- 
tions of the form 
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A,{x) = (l + xY''^ 



(23a) 
(23b) 



with 7i = 2.5 and 72 = 3. 
To test the general applicability of these expres- 
sions, we have run the simulations varying the dura- 
tion of the shock from r = 0.01 tdyn,h to r = Atdyn,hi 
and also varying the shock amplitude, lexp- AH of the 
simulation results scale with /gj.„, as expected. We 
find that for t < tdyn,h (on the lower end for the real 
disk shocking), equations ( p^ ) describe the results of 
the simulations very well. For the more prolongated 
shocks, the slope of the adiabatic correction becomes 



shallower. Nevertheless, the results can be fitted by 
the same functional form, with the exponents given 
in Table |. 

Note that for the large values of t, the overall effect 



of shocking decreases {lexp 



) and so does the 



accuracy of the fits. However, it is likely that the 
adiabatic corrections are less steep for t > tdyn,h- In 
the limit of "slow shocks" , the asymptotic Weinberg's 
result (nj) becomes valid, as expected. 

We have also checked whether the adiabatic cor- 
rections are dependent upon the structure of the clus- 
ter. We have run the self-consistent simulation for an- 
other King model, with the concentration parameter 
c = 1.5. The fitting formulae (cq. [E3[) slightly under- 
estimate the energy change and dispersion in the core 
and the extreme edge of the cluster, but they hold for 
most of the intermediate range. The statistical errors 
are small and do not affect the comparison. Still, it is 
possible that our fitting expressions might break for 
very concentrated clusters and thus should be taken 
with caution. 

Finally, Figure |lj shows the energy changes at the 
end of run D, the self-gravitating simulation of disk 
shocking. As in the impulsive case, the stellar energies 
are enhanced everywhere in the cluster. Our fit to 
the simulation results has no free parameters: it is 
a sum of (A_E) expected for the fixed potential (eqs. 
fig] and [g3[) and the self-similar potential change 
(eq. [£l|). Although not exact, the fits are reasonably 
good. Also, we find no additional contribution to the 
dispersion of energy due to potential fluctuations. 

4.3. Change of the total energy of the cluster 

One of the predictions of Weinberg's (1994) linear 
theory is that in the limit of large r, the change of 
the total energy of the cluster due to a tidal shock 
with a fixed amplitude is inversely proportional to its 
duration. Since in our simulations the amplitude of 
the shock scales as t~^, to conserve the integral of the 
perturbing force, we expect the total energy change 
to scale as AEtot oc t~^. 

We use the results of the simulations with various 
parameters r and lexp to test that prediction. Figure 
O shows that in the limit r ^ tdyn,h, the value of 
AEtot is constant. For the extended shocks, the en- 
ergy change follows a power-law. The fitting formula 
that combines both limits and provides a good match 



in the middle (Figure O) is 



AEtot = AEtot\r=0 1 



-3/2 



(24) 



dyn,h 



This assures the validity of the linear theory in pre- 
dicting the integral characteristics of the system. 

5. Test Case: Radial Shock 

The main effects of the evolution of the potential, 
an additional heating in the core and the decrease 
of the energy dispersion in the middle of the cluster, 
are new. To check that these effects are not restricted 
only to the one-dimensional disk shocking, we investi- 
gate another type of tidal perturbation, radial shock- 
ing, and use a completely independent numerical code 
to examine this phenomenon. The instant accelera- 
tion is proportional to the radius-vector of the star 
and is directed towards the center: 



dt " ' 



imp ■■ 



(25) 



The symmetry of the perturbation in a sph erical clus- 
ter allows u s to use Spitzer's shell method (^pitzer & 
Hart 197l|) to test the SCF code. We wiU see that 



the qualitative properties of our solutions survive the 
change in numerical methods. 

5.1. The Shell method 

The shell method, originated from an early work by 



Hcnon (1964), makes use of the spherical symmetry 
of the system. A shell i represents an ensemble of 
stars at the radius r^ with the energy Ei and angular 
momentum Ji. The shell moves in radial direction 
with the velocity Vr,i and can freely penetrate other 
shells. All shells are assumed to have the same mass, 
Ms ■ For shells arranged in order of increasing radius 
Ti , the equations of motion are 



df^ 



J? (» - y2)GMs 



(26) 



The factor Y2 takes into account self-gravity of the 
shell i. While the equations of motion in the spherical 
system enjoy the simplicity of Newton's theorem (the 
force on a shell is determined entirely by the number 
of enclosed shells), the calculation of the potential 
involves all shells: 



$,; 



GMJ 



Ns 



E 



GAL 



(27) 



where N^ is the total number of shells. Therefore, the 
energy of the shell, Ei = ii^ j/2 + Jf/2r^ — <I>j, depends 
not only on its own motion, but also on the position 
of the shells outside. 

This method was successfully applied by Spitzer 
and collaborators at Princeton to study the evolution 
of a spherical star cluster under the influence of the 
explicitly imposed perturbations imitating two-body 



encounters between stars (see Spitzer 1987 for a de- 
tailed reference). A small energy imbalance resulting 
from the perturbations was corrected for by adjusting 
the kinetic energy of the shells. 

Testing the method without external perturba- 
tions, we noticed that the total energy of the sys- 
tem was not conserved. The energy grew monotoni- 
cally with time with the rate proportional to the time 
step. The reason for the nonconservation of energy 
is the force discontinuity when two shells cross each 
other. Accordingly, we introduced a two-step correc- 
tion procedure. First, instead of solving the second 
order differential equation ( p6| ) as a pair of the first 
order equations, we integrated the equation once to 
obtain a first integral of motion, Cf. 



drj 



Jf {2i - l)GMs 
72 T 



Ci = const. (28) 



This equation was then solved numerically to advance 
shells in radial direction. The advantage of solving 
this equation is in continuity of the velocity v^i dur- 
ing shell crossing. The second correction takes care 
of readjusting d after each shell crossing. This pro- 
cedure assures conservation of the total energy of the 
cluster, Etot, which has to be a linear combination 
of the integrals of motion, Ci. As is easily checked, 
Etot = X] C'j/2, with the sum extending over all shells. 
Our test runs show exact conservation of Etot in the 
absence of perturbations. 

Another correction is necessary to solve equation 
(28). When a shell is close to its turnover points, 
Tmax or rmin, in & given time step it may jump over 
to the region of imaginary velocities Vr ■ The positions 
of such shells are integrated in two steps, first up to 
the turnover point, and then back with the reversed 
velocity. 

We used the same initial conditions as for the 
three-dimensional simulations, projected onto the ra- 
dial direction. The angular momentum of the shells 
was calculated from their initial positions and veloc- 
ities and then kept constant for the rest of the sim- 
ulation. We have checked that the initial conditions 



represent the assumed King model with correct den- 
sity profile and energy distribution. 

The nature of the correction procedure for shell 
crossing is intrinsically serial, due to the simultane- 
ous adjustment of the constants Ci for both shells 
involved in the crossing. The number of crossings 
in a given volume grows very fast with iV^, inflating 
the number of required operations. The largest run 
reasonable in computational time involves 10^ shells. 
This disadvantage is compensated for by the high pre- 
cision of the results. 

Figure n3 shows the energy changes for the impul- 
sive radial shocking (run E). This plot is qualitatively 
similar to the impulsive disk shocking (see Figure S). 
We are again able to fit the results with the sum of 
the analytical prediction for the fixed potential and 
the self-similar potential change: 



(A£;) 



2 

(Ai?^> = ^^H^r%L.(l + X..).(29b) 



2 r^-fc/Lp (-$,), (29a) 



This simple fit seems to hold for the radial shocking 
as well. However, once again the energy dispersion is 
decreased by the postshock evolution. The effect at 
the end of the simulation is even stronger than it was 
in the case of disk shocking. This negative feedback 
must be due to an even stronger anti-correlation of the 
effects of the fluctuating potential with the motion of 
stars at the time of shocking. 

5.2. Comparison with SCF simulation 

For comparison with the shell method, we have run 
the three-dimensional SCF simulation with the radial 
perturbation (eq. p5||). Run F involves the same 
initial conditions as the disk shocking case. 

The results of the 3D simulation and the ID shell 
method are very similar (cf . Figures 11^ and 11^) . The 



same fit (eq. 1 29 ) is used for both runs and describes 
(A£') adequately. The reduction of the dispersion is 
present in both simulations, constituting a convincing 
result. We have also checked that numerical errors do 
not contribute to our results. 

6. Summary 

We have explored tidal disk shocking on a star 
cluster using two independent types of self-consistent 
N-body simulations. As expected, the effects in the 



outer parts of the cluster are well described by the 
classic impulse approximation, but at the half-mass 



A detailed Fokker-Planck treatment of the clus- 



radius and in the inner parts, the effects arc more 



complex. 

One of the main goals of this paper was to investi- 
gate the conservation of adiabatic invariants of stellar 
orbits and to obtain convenient adiabatic corrections 
that could subsequently be used in the Fokker-Planck 
calculations. Solution of the Fokker-Planck equation 
is easier and much faster than full N-body simulations 
(see, for example, Spitzer 1987). The details of our 
F-P code are given in Gnedin, Lee, fc Ostriker (1998)|. 



ter evolution will be reported in Gnedin, Lee, fc Os- 
triker (1998)| . But we can note here what changes 



The Fokker-Planck code assumes a Gaussian distribu- 
tion of the energy changes of individual stars around 
the mean value, (AE'), in each energy bin. We have 
checked that the final distribution in our simulations 
is very nearly normal (Figure |7|), allowing us to use 
the two variables, (AE) and {AE'^), to describe the 
effects of shocking. Alternatively, one can use the full 
for malism of the linear t heory of Weinberg (1994a-c) 
and Murali fc Weinberg (1997a-c). 

First, let us address the mean energy change, 
{AE). It is given by equation (Pa|), where we find 
the adiabatic correction Ai{x) given approximately 
by (1 + x^)^^!"^ for shock durations comparable to 
the half-mass dynamical time, and by (1 -|- a;^)~^/^ 
for "slow shocks" . These results apply when the post- 
shock evolution of the potential is neglected. If the 
latter is included, the net effect is to produce an addi- 
tional heating due to the self-similar potential read- 
justment, AEpot = c/j^„p(— $i), where c is a fixed 
normalization constant (not a free parameter!). This 
second change is qualitatively similar to the method 
used by Gnedin fc Ostriker| (1997) for correcting the 
cluster energies to allow for adjustment to a new equi- 
librium, but different in detail in that the exact treat- 
ment allows for the exchange of energy between the 
stars, which occurs after the shock. 

The second order energy changes have similar adia- 
batic corrections (cf. Tabled), varying from (1-f x^)~^ 
to (1 -I- x^)"^/'' as the shock duration increases. In a 
self-consistent (time- varying) potential, the postshock 
evolution tends paradoxically to decrease the energy 
dispersion because of an anti-correlation between the 
initial energy change and the subsequent virialization. 

Finally, we found that the change in the total en- 
ergy of the cluster as a function of shock duration is 
in good agreement with the prediction of the linear 
theory. When the integral of the perturbing force is 
kept fixed, AEtot oc t"^ for r > tdyn,h- 



are expected when the postshock oscillations are al- 
lowed. We expect that incorporating the additional 
first order effect {{AE)) into the Fokker-Planck code 
will increase the rate of cluster evolution, leading to 
more rapid core collapse. Paradoxically, the second 
order effect ((A_E^)) of postshock oscillations appears 
to lead to a small reduction in the overall value of the 
relaxation and a consequent reduction in the rate of 
the relaxation-driven evolution. 

We would like to thank Lyman Spitzer Jr., Jeremy 
Goodman, Douglas Heggie, David Spergel and Scott 
Tremainc for discussions, Antonio Peimbert for sug- 
gesting an efficient implementation of the shell method, 
and the referee, Martin Weinberg, for a number of 
useful comments. We are especially indebted to 
Kathryn Johnston and Lars Hernquist for their help 
at the beginning on this project. Finally, we are grate- 
ful to Lars Hernquist for letting us use his original 
version of the SCF code. This project was supported 
in part by NSF grant AST 94-24416. 



A. Parallel implementation of the SCF code 

We have implemented the self-consistent field code 
on two different platforms, the SGI Power Challenge 
with shared memory architecture and the IBM SP2 
array with distributed memory. 

The benchmark tests on R8000 processors of the 
Power Challenge show that the computational time 
scales perfectly linearly with the number of particles, 
from N — 10* to 10^. The speed-up factor per pro- 
cessor is also close to unity, indicating that the code is 
ideally parallelizible. Figure |l^ illustrates the speed- 
up gain per time step per particle as a function of the 
number of CPU. The gain factor is about 3.8, 6.6, and 
11.5 for 4, 8, and 14 CPU, respectively. 

The performance of the SCF code on SGI Power 
Challenge can be compared to its implementation on 



a Connection Machine 5 reported by Hernquist et al 



1995) . The number of operations per particle de- 
pends on the number of coefficients retained in the 
expansion of the potential-density pair. For our test 
simulations, we have taken Umax = 6, Imax = 4. An 
independent measure of the performance is the com- 
putational time, t„, per time step per particle per 



10 



number of the expansion coefficients times the num- 
ber of CPU. [Hernquist et al. (1995) report i„ « 0.9 /iS 
for their nms (their Table 1). On the Power Chal- 
lenge, we get a little better speed, t„ = 0.65 fis. The 
average time per time step on a single processor in 
the run with 10^ particles was 120 s. 

We have also done a set of runs on SGI Origin 2000 
machine with RIOOOO processors. On a single CPU 
of Origin 2000, the code runs 2.9 times faster than 
on the corresponding R8GGG processor of the Power 
Challenge. However, the speed-up rate with number 
of CPU was a little lower. 

Finally, we have implemented the SCF code on 
IBM SP2 at the Maui High Performance Comput- 
ing Center using Message Passing Interface (MPI). 
Our preliminary tests show similar performance of the 
code on SP2 and the Power Challenge. 



B. The Position Velocity Correlation 

We calculate here the correlation function of the 
mean positions and velocities of stars of the same en- 
ergy, E. This function appears in the analytical ex- 
pressions for the energy dispersion due to tidal shocks 
(eqs. @ and H). 

Consider a spherically-symmetric stellar system. 
We seek to establish a relation of the form 



{r'v^)^^{r^)^{v')^{l+Xr.v), 



(Bl) 



-'0,0 

and therefore 



^2,0 

Io,o 



(r\') 



^2,2 
-^0,0 



1 + Xr,v 



^2,2 ^0,0 
-^2,0 -^0,2 



(B3) 



(B4) 



The integration over the velocity space is straightfor- 
ward with the aid of the delta-function. Also, the dis- 
tribution fmiction f{E) drops out of integration since 



the energy is constant. Cancelling some constants, we 
can redefine the interesting integrals as 



where the correlation function Xr,v{E) can depend 
only on the energy. Introducing the integrals 



= (47r)^ / / f{E)d[E- — - $(r) ) w^+"rf?;r^+™dr, 

(B2) 
where S{x) is the delta- function restricting the inte- 
gration to the stars of energy E, we can write the 
phase space averaging as 



^n.An, — 



J2+m 



[E - <^{r)]-^ dr, 



(B5) 



where r„i [E) is the maximum radius a star of energy 
E can reach, determined by the following condition: 
E = $(r,„). After a simple dimensionalization of the 
integrals, we finally define 



1 



<i>(xr„ 



$0 



$(r„ 



$0 



dx, (B6) 



where $o is an arbitrary constant. Note that the ex- 
pression for 1 -|- Xr.v has an equal number of the n and 
m indices in the numerator and the denominator, so 
that we can cancel from the integration factors such 

l + r. 



as r^™ and {2E - 2<i>o)~2~. Equation (B6) can now 



be solved numerically for any given potential <l>(r). It 
is useful to choose the value of the constant $o at the 
center of the cluster. 

Note, that if the potential obeys a simple power- 
law, <I>(r) = $0 + cr", the integral (B6) can be done 
analytically. The correlation function assumes a par- 
ticularly simple form. 



Xr,v 



10 + 3a' 



(B7) 



depending only on the power-law index, a. In general, 
Xr,v is always negative for self-gravitating systems. It 
varies from Xr,v = —0.25 for the harmonic potential 
to Xr,v ~ —0.57 for the Keplerian potential, and lies 
somewhere inbetween for other systems. 

We have calculated the correlation function for a 
number of King models with various concentrations, 
c, and fitted the results as a function of the normalized 
stellar energy e = \E/Eund\- 



1-l-X,.,^ = 0.75-0.1 (1+c) (l-e)-O.l ■ 



(1.8 -c)ci/" 
l + (e/0.15-l)2' 
(B8) 
The fit is accurate to 2% for the range of concentra- 
tion parameters c = 0.6 to 2.7. 

By doing the angular integrals, it is straightforward 
to show that for spherical systems with isotropic ve- 
locity distribution any combination of one-dimensional 
coordinates and velocities has the same correlation 
function. For example. 



{z^vDj^ = (Z^)^ {vl)j^ (1 + Xr,v)- 



(B9) 
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Table 1 
Parameters of N-body simulations 



Run 


Shocking force 


Direction 


Cluster potential 


A 


impulsive 


vertical 


fixed 


B 


impulsive 


vertical 


self-consistent 


C 


time- varying 


vertical 


fixed 


D 


time- varying 


vertical 


self-consistent 


E (shells) 


impulsive 


radial 


self-consistent 


F 


impulsive 


radial 


self-consistent 



Table 2 
Exponents of the adiabatic corrections 



Regime 



71 



72 



'^ ^ '^dyn,h 
'^ — ■^ ^dynji 



2.5 

2 

1.5 



3 

2.25 

1.75 



REFERENCES 

Binney, J., & Tremaine, S. 1987, Galactic Dynamics 
(Princeton: Princeton Univ. Press) 

Gnedin, O. Y., & Ostriker, J. P. 1997, ApJ, 474, 223 

Gnedin, O. Y., Le e H. M., fc Ostrike r, J. P. 1998, 
ApJ, submitted; |astro-ph/9806245 



Henon, M. 1964, Ann. d'ap., 27, 83 

Hernquist, L., & Barnes, J. E. 1990, ApJ, 349, 562 

Hernquist, L., & Ostriker, J. P. 1992, ApJ, 386, 375 

Hernquist, L., Sigurdsson, S., & Bryan, G. L. 1995, 
ApJ, 446, 717 

Johnston, K. V., Hernquist, L., & Weinberg, M. 1998, 
ApJ, submitted 

King, I. R. 1966, AJ, 71, 64 

Kruskal, M. 1962, J. of Math. Phys., 3, 806 

Kundic, T., & Ostriker, J. P. 1995, ApJ, 438, 702 



Landau, L. D. & Lifshitz, E. M. 1988, Mechanics, 3rd 
ed. (Oxford: Pcrgamon) 

Murah, C., & Weinberg, M. D. 1997a, MNRAS, 288, 
749 

Murah, C., & Weinberg, M. D. 1997b, MNRAS, 288, 
767 

Murah, C., & Weinberg, M. D. 1997c, MNRAS, 291, 
717 

Ostriker, J. P., Spitzer, L. Jr., & Chevalier, R. A. 
1972, ApJ, 176, L51 

Spitzer, L. Jr. 1987, Dynamical Evolution of Globular 
Clusters (Princeton: Princeton Univ. Press) 

Spitzer, L. Jr., & Hart, M. H. 1971, ApJ, 164, 399 

Tremaine, S. 1996, private communication 

Weinberg, M. D. 1993, ApJ, 410, 543 

Weinberg, M. D. 1994a, AJ, 108, 1398 

Weinberg, M. D. 1994b, AJ, 108, 1403 



12 



Weinberg, M. D. 1994c, AJ, 108, 1414 



* 






r 






''\ \ 


//■■■■' 
//■■ 



Fig. 1 . — Schematic picture of the readjustment of the 
cluster potential. Dots mark the initial King model, 
solid line represents new equilibrium configuration. 
(Radial axis is reflected for the illustration purpose.) 



This 2-colurnn preprint was prepared with the AAS lATJiiX 
macros v4.0. 
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Fig. 2. — Comparison of the Spitzer [A{x) = e~^^ ] 
and Weinberg [A{x) = (1 + a;^)~^/^] adiabatic correc- 
tions, i. e., the reduction of the energy change due 
to tidal shocks because of conservation of adiabatic 
invariants of steUar orbits, versus the adiabatic pa- 
rameter, X [eq. (||)]. For detailed self-consistent sim- 
ulations, we find that adiabatic corrections show a 
power-law behavior intermediate between these two 
curves (cf. eq. 12^ and Table 0) . 










Fig. 3. — Integral characteristics of the cluster af- 
ter the impulsive disk shocking (run B). Upper panel: 
the change of the total energy of the cluster (solid 
line) versus time. Dots and dashes show partial con- 
tributions of the kinetic and potential energies, re- 
spectively. Energy is normalized to the binding en- 
ergy of the cluster (i.e., the absolute value of the cen- 
tral potential), and time is in units of the dynamical 
time at the half-mass radius. Dotted horizontal lines 
mark the expected values for the kinetic and poten- 
tial energies after the virialization: ATvir — —AEtot, 
AWvir — 2AEtot- Lower panel: virial ratios —2T/W 
(solid line) and —2T/V (dashed line) versus time, 
where V is the Virial of the system (eq. [0|). On this 
plot both lines coincide. Horizontal dashes mark the 
1% deviation from the equilibrium, which we consider 
as the "acceptable" range. Virial oscillations last for 
about ten dynamical times. 
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Fig. 4. — Energy change and dispersion versus ini- 
tial energy of stars for the impulsive shocking in the 
fixed potential. Energy along the abscissa axis is nor- 
malized to the binding energy of the cluster. Data 
points show the energy changes at the end of the sim- 
ulation. The solid lines are the analytical predictions 
(eq. III3I). The center of the cluster is on the left; the 
outer parts are on the right. Vertical arrows indicate 
the average energy corresponding to the core radius 
and the half-mass radius of the cluster. The initial 
tidal radius is near the right edge of the plot. 




Fig. 5. — Energy change for the impulsive shocking, 
including self-gravity. Open boxes show the results 
immediately following the shock, solid triangle - at 
the end of the simulation. Dashed lines show the an- 
alytical estimates for the fixed potential (eq. ||l8||). 
{A.E) is enhanced everywhere in the new equilibrium 
state. The solid line is our fit (eq. [Q). Dispersion is 
decreased in the middle of the cluster; the apparent 
flattening in the core is due to numerical relaxation 
(see Text). The dashes on the upper panel show the 
estimate of the energy dispersion at the end of the 
simulation due to two-body relaxation. 
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Fig. 6. — Density profile of the cluster at the end 
of simulation B (histogram). Solid line is the initial 
King model. An apparent peak in the leftmost bin is 
an artifact of binning. 



Fig. 7. — Distribution of stellar energies in one bin at 
the half-mass radius for run B. Thin histogram shows 
the relative energy change immediately following the 
shock, thick histogram shows the final distribution. 
Both energy changes are normalized to the dispersion 
at the end of the simulation. Dashed and dotted lines 
are the Gaussian fits to the histograms. Note (1) the 
good fit of the final distribution to the Gaussian form, 
and (2) the narrowing (cooling) of the distribution in 
the final state. 
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Fig. 8. — Energy change due to potential evolution 
versus initial radial velocity in the bin at the half- 
mass radius for run B. Both quantities are in the code 
units. The postshock energy gain correlates with the 
motion of stars at the moment of shocking. 



Fig. 9. — Energy change due to potential evolution 
versus impulsive energy kick, A£^sh, in the bin at the 
half-mass radius for run B. Here (AEsh) is the aver- 
age energy change in the bin. All energies are in the 
code units. The anti-correlation of AEi — AEgh and 
Ai<^sh leads to the reduction of the dispersion after 
the virialization. 
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dyn.h 



Fig. 10. — Same as Figure 0, but for the time- varying 
perturbation (run D). Vertical dotted line across the 
plot at i = 4:tdyn,h marks the maximum amplitude of 
the shock. Lower panel shows virial ratios: —2T/W 
(solid line) and —2T/V (dashed line). 





Fig. 12. — Same as Figure |T^, but for the case of 
the self-gravitating potential (run D). The solid line 
in the upper panel is the analytic expectation for the 
fixed potential, the solid line in the lower panel is our 
simple fit including the potential energy change (eq. 
pl[). Dashes on the upper panel show the estimate 
of the energy dispersion at the end of the simulation 
due to two-body relaxation. 



Fig. 11. — Same as Figure^, but for the time- varying 
perturbation in the fixed potential (run C). Dots show 
the analytical expectation with the Spitzer adiabatic 
corrections (eq. [|l^), dashes correspond to the Wein- 
berg corrections (eq. JMJ ) , soHd lines are our fits (eq. 
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Fig. 13. — Change of the total energy of the clus- 
ter due to time-varying shocking as a function of its 
duration for the runs with various shock amplitudes. 
Dotted lines are our fits (eq. [p4J). 
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Fig. 15. — Same as Figure |5|, but for the radial impul- 
sive shocking (run F). The fits are given by equation 
(P9|). Stars show the numerical errors at the end of 
the simulation as estimated in the isolated run with- 
out perturbation. The errors are negligible except for 
the dispersion in the core of the cluster. 



Fig. 14. — Energy change for the impulsive radial 
shocking using the shell method (run E). Dashes are 
the analytical estimates for the fixed potential, solid 
line in the lower panel shows the fit including the 
potential energy gain (eq. [p9[). 



19 




5 10 

number of CPU 



Fig. 16. — Scaling of the computational speed-up per 
particle per time step with the number of CPU on 
SGI Power Challenge. The points represent test runs 
with various numbers of particles, from N =- 10^ to 
10^. The diagonal dashed line shows the ideal scaling 
law. 
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